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ABSTRACT 


A three-dlmeneioiial, time-dependent mathematical model to 
represent Mobile Bay was developed. The objective of this study 
was to develop computer programs which would numerically solve the 
appropriate conservation equations for predicting bay and estuary 
flow fields. The model will be most useful for analyzing the dis- 
persion of sea water Into fresh water and the transport of sediment. 
Also, the model serves as a useful tool for relating field and phy- 
sical model data. The unique feature of this model Is that It 
correctly accounts for momentum transfer In the governing flows; 
thereby, making It far more realistic than any previously devised. 
NASA's KRTS and Sky lab programs resulted In high quality photographs 
of Moi>lle Bay. U.S, Army Corps of Engineers have also studied this 
bay e7;tenslvely. All these data have been reviewed for comparison 
to this mathematical model. 
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INTRODUCTION 


This report describes the development ot’ a three-dimensional, 
time-dependent flow field model to represent Mobile Bay. This model 
is a mathematical representation of the appropriate conservation 
equations which has been evaluated with a digital computer. Since 
all important physical phenomena cannot be mathematically modeled 
in a definitive fashion, the errors introduced by approximations 
made to make the model tractable must be evaluated. Therefore, the 
model was prepared specifically to represent Mobile Bay. 

Mobile Bay is not only a commercially Important but also a 
technically Interesting bay to study. The bay is connected by a 
tidal inlet to the Gulf of Mexico, thereby forming a region in which 
sea water is measurably diluted with land water drainage. Such a 
region is called an estuary. The drainage occurs by several large 
rivers discharging both water and sediment into the bay. A ship 
channel extends from end-to-end of the bay; the channel is dredged 
and maintained by the Corps of Engineers. The bay is shallow and 
therefore greatly affected by prevailing winds. 

All estuaries are so complex that any mathematical model which 
may currently be implemented on a computer will contain certain 
simplifications to make a solution of the model tractable. These 
simplifications are so severe that once the model is developed it 
must be tested to prove that it is a valid representation of the 
estuary. Ideally, one would like to have sufficient experimental 
data from direct field measurements to test the model. Usually, 
sufficient data are not available. Even when they are, initial 
attempts at modeling will Invariably show some discrepancy, then 
the model will have to be "improved". 

Methods of making such improvements are not, usually, obvious, 
because the very complexity of estuaries defies establishing cause 
and effect phenomena. Results from physical model studies contri- 
bute much more to this stage of mathematical model development. 

This is due to the fact that physical model experiments can be 
designed to emphasize specific study of select phenomena, i.e. com- 
plicating factors can simply be eliminated in these studies. Use 
of only "real world" field data cannot be so utilized because our 
environment cannot be so controlled. Of course, physical model data 
is also often used without collaborating field or mathematical model 
data. 


These general comments concerning the interactions of field, 
mathematical model, and physical model data precisely describe our 
current state-of-knowledge of Mobile Bay. Several attempts have 
been made to measure the flow properties within the bay; the most 
complete and recent study was performed by the Corps of Engineers 
to provide verification data for a physical model study which is 
being conducted, also by the Corps at Vicksburg, Mississippi. The 
purpose of the physical model experlmentF is to study the Impact of 
opening another ship channel within Mobile Bay. Photographic data 
from NASA are also available. 
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ThlB report summarizes the Corps of Engineers data pertinent 
to Mobile Bay, describes the mathematical model which has been 
developed to describe Mobile Bay, and conipares data from all three 
sources. 


SUMMARY OF FIELD lATA 


Results from several field studies of Mobile Bay have been 
reported (1, 2). Although these reports present many measurements 
descriptive of the flow-field, there is not sufficient information 
at any given instant of time to constitute adequate boundary condi-. 
tions and check-points to be compared to a mathematical model. An 
adequate experiment was performed by the Mobile District Corps of 
Engineers for verifying the physical model study which will be 
described in the next section of this report. 

These verification measurements (3) consisted of flow and tide 
level data through the four major river inlets to the north and the 
two major passes in the south for a complete tidal cycle. At the 
passes, vertical salinity profile data were also collected, Cross- 
8‘“ctional areas for the inlet and outlets were determined. Finally, 
tide, velocity, and salinity vertical-profile data were also taken 
at select stations within the bay. These regions are shown in 
Figure 1. 

There are other rivers which were not measured during this 
study, presumably because there is little flow through them. How- 
ever, what was measured corresponds exactly with the streams which 
are included in the physical model. Also, no data on winds or on 
flow into or out of the northern marshes are given. There is an «• 
airport weather station nearby from which some wind data may be 
obtained. 

The four river inlets were taken to be where the highway crosses 
the north end of the bay. These inlets and cross-sectional areas 
were the; Mobile, 2,120 sq, m (22,810 sq. ft,); Tensaw, 2,850 sq, 
m (30,660 sq. ft.); Apalachee, 2,013 sq, m (21,660 sq, ft.); and 
Blakeley, 2,499 sq. m (26,880 sq. ft). These rivers are shown in 
Figure 2, Notice that part of the flow from the upper Mobile River 
passes into the bay through the Tensaw inlet, and that all of the 
flow from the upper Tensaw enters through the Tensaw, Apalachee, 
and Blakeley. 

The mouth of Mobile Bay is 30,265 sq. m (325,600 sq, ft.) in 
cross-section. The pass connecting Mobile Bay to the Mississippi 
sound between Cedar Point and Dauphin Island is 6,144 sq. m (66,100 
sq. ft.) in cross-section. All cross-sections are the area between 
mean-sea level and the bottom. The bottom location changes signi- 
ficantly with time, probably because of sediment deposition. 

The previously mentioned data are available for further study. 
Only two points need further commment at this time. Lov? tide at 
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Mobile Point occurs some 4 hours before low tide at the Mobile 
River inlet; high tide about 2.5 hours before. This demonstrates 
that there tire significant tidal phenomena within the bay. 

The second point is that the high tide causes flow reversal 
in the river Inlets at the north end of the bay. This implies that 
one would have to go up-stream before the rivers would be discharg- 
ing, l.e. the bay, marshes, and lower river channels are filling 
with water on a rising tide. This fact complicates the mathematical 
specification of boundary conditions at the northern end of the bay. 
These marsh regions ara clearly discernible in the ERTS IR photograph 
shown in Figure 3. The dark regions in the river valleys north of 
the bay are the marshes. Notice that these marshes do not extend 
down either side of the bay. 

F^^ure 4 is a companion picture to Figure 3, but it was taken 
at visual wavelength. Although quantative data cannot be obtained 
from such photographs, the long distances over with the river plumes 
retain their shape and the dead water region in Bon Secour are 
clearly evident. Similar data in true color from Skylab photographs 
are also available. 


PHYSICAL MODEL lATA 


In order to estimate possible effects caused by building a new 
ship channel, to be named Theodore Channel, the Corps of Engineers 
have performed physical model studies of Mobile Bay. First a set 
of base data for the Bay as it now exists was collected, then the 
new channel and several possible new islands which serve as spoil- 
age banks were modeled and studied. The base data are of interest 
to this study; such data were made available to these investigators 
by the Mobile District Office of the Corps (4). 

The region of the Bay which was physically modeled is shown 
in Figure 1. Salinities and velocities at the numbered surface 
locations during a tide cycle are reported for the "surface" and 
"bottom" elevations. The Mobile and Tensaw Rivers flows were 
metered and reported, presumably well up-stream of the northern 
reaches of the Bay. Note this does not correspond to the location 
where the field data are reported. The river volumetric flow rates 
were reported as (cubic feet/second) of prototype flow, and tide 
levels in (ft) of prototype surface elevation. The scaling laws 
used for such unit conversions are not presently known to these 
investigators. Note also should be made that "turbulence strips" 
were used in these model studies. Therefore, the physical model 
data will be considered as its representation of the actual Bay, 
which is how the modelers intended it should be used. 

Tide data at the north, south, and middle points of the Bay 
are shown in Figure 5. Note that high tide occurs about 3 hours 
later at the state docks location than at Dauphin Island. This 
difference between physical model and field data is thought to be 
significant. 
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Figure 3. IR Photograph of Mobile Kay from NASA's KRTS Program 
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Figure 4. Black and White Photograph of Mobile Bay from 
NASA's ERTS Program 
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All velocltlee are reported as north or south velocity com- 
ponents only. Data are available at many more model stations than 
fleldl stations. Interior stations at B-12 and B-32 are the only 
two field stations. 


THE MATHEMATICAL MODEL 


The philosophy of developing and using a mathematical model 
Is that a certain set of equations represents the phenomena of 
Interest In the region of interest and that these equations may be 
solved for a particular set of auxiliary conditions. The solution 
would represent the behavior that one wishes to predict, and the 
auxiliary conditions and/or specified parameters In the equations 
would represent alternatives that one wishes to either understand 
or control. 

Any and all errors which might be In the model would manifest 
themselves as Inaccuracies In the flow description, hence, calcu- 
lated results should be verified before a model Is used. All known 
Information concerning the flow region of Interest should be con- 
sidered to establish this verification. It Is Impossible to either 
completely measure or model a large, complex flow- field; therefore, 
model development and verification proceed in a "hand- over- fist" 
fashion. That Is, as more is learned of the flow better predic- 
tions can be made, and as better predictions are made the more 
important regions of the flow can be Identified for further field 
study. 

Any model study Is bounded by the desire for completeness and 
the necessity for econonqr; therefore, a modeling scheme which 
attemps to optimize this balance was developed. Since flow within 
the ship channel Is believed to exert a dominating Influence on 
the flow of sea water Into the bay and on the sediment deposition 
trends within the bay, a three-dimensional, time dependent flow 
field model was deemed essential to describe the type flow found 
In Mobile Bay. 

Model Development 

All solutions to fluid mechanics problems are based on conser- 
vation principles. These principles, usually expressed as a con- 
servation of mass, momentum and energy, are presented In most stan- 
dard texts of fluid mechanics for single-component flows. The 
phenomena of interest Is that Involving the relative motion and 
the mixing of fresh and saline water. Sediment transport Is also 
of Interest, but its study will be postponed to a future time. 

Because the fluid In the mixing region may have two components 
some care must be exercised to Insure that the most useful form of 
the pertinent conservation laws is determined. Also, It is fruit- 
ful to consider and compare experimental results and analytical 
methods which arise from studies of heated water discharges. The 
relationship between the mass density of fresh water and sea water 
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and various mixtures of the same may be represented as a linear 
function of the amount of salt in the mixture, but the number of 
moles (or particles) in a unit volume is a constant. This means 
that any diffusion which occurs must behave as a replacement process, 
i.e. a molecule of ’’salt-water" replaces a molecule of fresh water. 
Furthermore, by introducing a fictitious molecular weight for salt 
water, the fluid may be described as a binary mixture of fresh 
water molecules and salt-water molecules. This behaviour of the 
number density allows one to specify the condition for incompres- 
sible flow by stating that the substantial derivative of the number 
density of fluid particles with respect to a number average velocl.py 
is zero, or that the divergence of the number average velocity is 
zero for an isothermal fluid. A similar relationship for the effect 
of modest temperature changes on the density of water exists. Now 
all the conservation laws may be stated. 

Very generally, the conservation of mass equation states that 
the amount of mass accumulated in a given volume is eqx}al to the 
difference between the mass convected out of the volume and the mass 
convected in. For multicomponent fluids, the net flux of mass into 
the volume by diffusion must also be considered. 

The conservation of momentum is derived from Newton's second 
law of motion. Expressed for fluids, this states that the rate of 
change of momentum of the fluid in a given volume is equal to the 
summation of the vector forces acting on the fluid. Pertinent 
forces may consist of pressure forces, shearing stresses, and grav- 
ity forces. An apparent force known as the Coriolis force may also 
be included when the equation is expressed in^a coordinate system 
which is moving relative to an Inertially fixed system. Unfortun- 
ately, the momentum and energy equations are cumbersome when written 
in terms of the number average velocity; therefore, both the number 
average velocity and the mass average velocity will be used in 
formulating the conservation laws. 

The conservation of energy equation results from applying the 
first law of thermodynamics to the moving fluid; it states that the 
rate of Increase in total energy is equal to the sum of the rate 
of work done on the fluid and the rate of heat added from external 
sources. The total energy of a given volume of fluid consists of 
three types: internal energy which may be expressed as a function 

of the temperature, kinetic energy due to velocity of the fluid, 
and potential energy which is a function of elevation of the parti- 
cles of fluid. The rate of work done to the fluid results from 
pressure forces, gravity forces, and viscous and turbulent shearing 
forces. 

The energy equation reduces to a trivial form for incompressible, 
isothermal flows When velocities are such that there is little heat 
generated by viscous dissipation. Thus, the conservation of mass 
and momentum equations produce an independent and complete set. 

In addition to the general conservation equations, an equation 
of state is required to express density as a function of composition 
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and temperature. For liquids, density and temperature are usually 
Independent o£ pressure, l.e. Incompressible. 

The conservation laws and equations are summarized In Table 1. 
Note that Interpretation of the transport coefficients as either 
molecular or eddy values allows the stated equations to be used 
for both laminar and turbulent flows. The complete conservation 
equations just described represent a complex set of second order, 
nonlinear partial differential equations. Even so, these equations 
do not precisely describe certain types of flow. This Inexactness 
Is due to lack of properly specifying specific terms of these 
equations as certainly the basic laws must be observed. It Is an 
easy matter to Indicate turbulent transport with a coefficient; 
specifying the local value of such a coefficient is another matter 
entirely. Furthermore, the microscopic interaction of the fluid 
particles and solid material (l.e. sediment) is not completely 
known, and is not, therefore, represented by the equations in 
Table 1. 

The actual solution of these equations will use an approxi- 
mation presented by Frank-Kamenetskii (Ref 5) and utilized by 
Daly and Pracht (Ref 6) which states that negligible differences 
result from using the number average velocity instead of the mass 
average velocity in the conservation of momentum and energy equa- 
tions. This then yields a set of equations with only one velocity. 

Another approximation suggested by Yih (Ref. 7) is quite use- 
ful in simplifying the conservation equations, Yih recognizes that 
the coefficient of expansion dp/3T is veiry small compared with unity 
for most liquids and, through a perturbation analysis, shows that 
V • U s: 0 when the range of temperatures are small relative to the 
mean absolute temperature. The implication of this assumption is 
that the volume dilation due to thermal expansion is negligible. 

This approximation is analogous to that involving the use of 
number-average velocity to simplify the solution of the two-component 
system of equations. 

Before the equations are solved, one must choose a suitable 
coordinate system. A Cartesian system was chosen. 

The resulting equations, including the approximations described 
above, which apply directly to non-isotheinnal, two-component fluid 
motions for the case of interest are may be solved by the method 
reported herein. The transport co-efficients may differ in each 
of the spatial directions. It should also be noted that the radi- 
ative heat flux vector qr of the energy equation need only be con- 
sidered at the surface of the fluid; hence, it could be justly 
considered as a boundary condition. Furthermore, the two most common 
situations arise when the fluid is one-component (S = 0) or is 
isothermal (T = constant); seldom must both variations be considered 
8 imul tane ous ly . 

Although the species continuity equation of Table 1 is written 
for the salt concentration S, it is equally applicable for the 
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Table 1. General Conservation Laws for Modeling 
Bays and Estuaries 


Overall Continuity Equation. 

+ U . 7 p + p(7*U) “ ” V) 

ot 

Approximate Equation o£ State* 


o 4. (Ifil 

Pf ^ Us) 


+ (S CT - T ) 


Diffusive Behavior. 

Ill (If ^"-4 


p7*(U - V) = 


Compressive Behavior. 


p(7*U) = 

Momentum Equation. 
9(pyi 4 . v»(7pV) = 

Energy Equation. 


(g + U.Ttj,, 


TP + T.p.7V - pg -• p(2fbcV) 


(1) 


( 2 ) 


(3) 


(4) 


(5) 


+ 7*(pEV) + 7*q + T'Ct - IP)*V - 2g^*j^ = 0 

Relationship Between Mass-Average and Number-Average Velocity. 


( 6 ) 


U - V = 

Relationship Between Mass Flux Vector and Molar Flux Vector. 

Jg « Kg[l + (S/p)(M^-Mg)/Mg] 


(8) 


conservation of any chemical or biological species with the Inclu* 
slon of a rate of generation (or degeneration) term. These equa- 
tions could then be applied to the dispersion of any pollutant 
discharged Into a .river or bay. Should the species have a negligible 
effect upon the density as Is often the case, then the species 
continuity equation could be decoupled from the other conservation 
equations. For this case, the flow field could be computed and 
the resulting velocities used to compute the advectlon, diffusion 
and generation of the species. For cases where the density Is 
affected by the concentration of the species, this decoupling Is 
not applicable. The species continuity equation should then be 
solved simultaneously with the other conservation equations. 

To model Mobile Bay, the flow Is considered Isothermal and 
to be a binary mixture of salt and water; therefore, the equations 
given In Table 2 apply. Furthermore, the viscous terms In the 
z-momentum equation were also neglete<‘(i, since an order of magnitude 
analysis showed them to be two orders lower than the pressure and 
gravity terns In this equation. 

Solution Technique 

The set of equations as presented In Table 2 Is sufficient 
to describe all the dependent variables for a given set of eddy 
coefficients, but the equations are In a form for which there Is 
no known analytical solution. These equations will be approximated 
by a set of finite-difference equations and solved by using a time- 
dependent technique. Before writing the equations in finite-dif- 
ference fetra, a stretching transformation that will create a more 
efficient use of grid points is suggested. For example, by letting 



where Cj and Cg are arbitrary constants, the equations of Table 2 
can be expressed as functions of Y instead of y. The implications 
of this transformation can be more easily understood by rewriting 
It ds 

Y = Ca tan"^ (qy) (16) 

When y ® 0, then Y = 0. Grid planes located at even increments of 
AY will produce ever-increasing Increments of Ay. This tranforroa- 
tlon will increase resolution where It Is most needed and make 
specifying boundary conditions less critical. 

The philosophy of using this numerical procedure is to begin 
with a set of Initial values of the dependent variables at all grid 
Intersections and allow the flow field to adjust asymtotlcally 
until all the conservation equations and boundary conditions are 
satisfied. The final result will be a steady-state flow field. 
Naturally, the better the estimate for the initial data, the more 
rapid the convergence, but even a crude guess, one conveniently 
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Table 2. Specific Conservation Laws for the Proposed 
Mobile Bay Model 


Species Continuity Equation. 


DS 

Dt 


n c?S . 




n 


Equation of State. 



Compressive Behavior, 

7 . u.|!i. +|!iy +|H. 
ox 3y oz 


Momentum Equation, 


D(pU» ) _ , 3__ rHir wBUY 

Dt 3x + dx 


4- ^ 4-A 


ay 


ay 


+ 2pn,Uy 

D(p-Uy ) _ _ ^ J. ^ 


Dt 


ay ^ ax 
2pn*u* 


ax 


\ J 


dz 


az 






Dt 


3z 


P8* 


( 9) 

( 10 ) 

( 11 ) 

( 12 ) 

(13) 

(14) 
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specified, is sufficient. The boundary conditions will then be 
replaced with conditions which depend on time, and the asymptotic 
solution will become an initial condition. An unsteady solution 
will result from further calculations. An efficient method for 
specifying initial conditions when making a series of calculations 
with similar geometry is to store the results of one run on tape 
and recall these data as initial conditions for the next run. 
Experience has shown that this can reduce computation times for 
subsequent runs by a factor of four. 

The numerical technique for calculating the dependent variables 
actually involves a sequence of calculations for each time step. 

The first step of the sequence involves computing new values of 
Ux > Uy« and S (and T, when needed) at each grid intersection by 
stepping forward in time by At. This is best explained by consid- 
ering the following truncated Taylor series; 

Ux(n+l,l,j,k) = U*{n,i,j,k} +|^fn,i,j,k} At (17) 

The braces indicate functionality, n is the index of the time step, 
and i,j,k indicate specific x,Y, and z grid planes, respectively. 

A finite difference approximation to the transformed right-hand 
side of the x-momentum equation of Table 2 can be substituted for 
9Ux/9t. Similar expressions can be written for Uy and S or T. 

The maximum possible time step for this type of numerical 
scheme is approximately given by the Courant-Friedrlcks-Lewy (CFL) 
condition, but in actuality this maximtun time step is seldom 
obtained. There are various procedures for Increasing the stable 
time step of a truncated Taylor's series such as the one shown, 
but they all pay a penalty in computation time per step and/or 
computer core storage. A method which has proven quite useful to 
these authors is best explained by considering another truncated 
Taylor series: 

+ U8) 

The braces indicate functionality, n is the index of the time step, 
and ifj»k indicate specific x,Y, and z grid planes, respectively. 
When (S"Ux)/Ot") is approximated by a forward difference as a 
function of (9Ux)/(9t), one has 

UxCn+l,i,J,k} - Ux{n,i,J,k} ,-)• [|^{n,l,j,i} + 

|^{n+l,ij,k}]|^ (19) 
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The (5H,y(at) {n,i,J,k3 can eaolly be computed from a finite- 

difference approximation to Equation 12 , Upwind differencing was 

used for the convective terms. For the case of positive and 

negative U/ and components of velocity at grid point i,j«k, 

the finite difference approximation of Equation 12 becomes 
* 

|^{n,i,j,k} » -'[|;^(Cu,Mn,i+l,j,k3 - U,Mn,l,j,k3] 
pln,i,j ,'k) t!P{n,l+l, j ,i] - PCn,i-l,j,k}3^ 

+ ^Ju*Cn,i,j,k} iv{n,l,j,k} - U, Cn,l,J-l,k3 [n,i, j-l,k}) 

4-^ (Ux{n,l,j,k3 Uxtn,i,j,k3 - U« {n,i, j ,k-l3 U* {n,i, j ,k-l3 )] 

+ Ch ■ 2Ux{n,i,J,k3 + U* {n,i-l, j ,k3) 

+ (UxCn,l,j+l,k3 - 2Ux{n,l,j,k3 + Ux{n,l,J-l,k3) 

+ (Ux{n»l.J+l,k3 - Ux{n,i,j-l,k3)] 

+ (Ux{n,l,J,k+l3 - 2Ux{n,l,j,k3 + U* (n,l,J,k-l3)] (20) 

where 

V u) = f {j} 

and 

V" tj} = 0 (j) 

Computing dUx /dt{n+l,i, j,k3 is more difficult because depen- 
dent variables at time n+1 are required but are not known. For 
this reason, solutions of Equation 19 usually fall into two classes: 
the implicit sahemes which Involve matrix inversion, and two-step 
schemes which estimate provisional values at time step n+1. From 
a practical stand-point, the primary reason for Including (SPUx)/ 

) term of Equatio'' 18 is to dampen numerical oscillations of 
(8Ux)/(3t) by averat+nr' gradients at time n and n+1, as shown in 
Equation 19. This effect may also be effectively accomplished by 
averaging gradients at times n and n-1, which would give 

Ux(n+l,i,j,k3 = Uxtn,i,j,k3 + [|^Cn,i,j,k3 + |^ [n-l,i, j ,k3 ,]|^ (21) 
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Equation 21 was, therefore an efficient method for ccrmimtlng 
U,{n+l,i, j,k3. Since (3U, )/(5t) {n-l,i,J,k} was already available, 
less manipulation per time step was required than for other commonly 
used procedures. Equation 21 provided stability comparable <o that 
of a two-step technique* Both techniques have time step limitations 
as determined by the Courant-Friedericks-Lewy conditions., 

It is conceded that using Equation 21 instead of Equation 19 
possibly sacrifices some accuracy in describing unsteady behavior, 
but if a steady-%tate solution is the desideratum, then lagging 
(3u)/(dt) by half a time step is Inconsequential because all partials 
with respect to time should asymptotically approach zero. 

Once new values of U, and Uy have been computed at each grid 
Intersection Uj may be calculated from a finite-difference represen- 
tation of equation 11, Spatial marching can be used in the vertical 
directioti, beginning at the bottom, where U, » 0. 

The free surface height h at each horizontal grid intersection 
may then be calculated. The incompressibility of the fluid requires 
that the net masa of fluid convected into a volume of fixed horizon- 
tal dimensions near the surface has to be accompanied by a comparable 
adjustment to the surface height of the volvaie. Since U, , Uy and 
at time n+1 have already been computed in the previous sequence, 
sufficient data are available to calculate (Sh)/(5t) {n+l,i,j}. The 
new surface height can be computed from: 

h{n+l,i,j} = hCn,l,j} + [||Cn,l,J} + |~[n+l,l,J}] (22) 

The gradient Oh)/(9t)(n,i, j } will be known from the previous time 
step. Note that h is related to the size of the grid cells which 
contain the surface; therefore, these cells change volume during 
the course of ‘the calculations. 

The size of the surface elements may be accounted for by iden- 
tifying (h - Zb#,jc-i) = ® then writing the divergence and momen- 
tum equations as follows: 




zmax-1 ^ Aa _ q 
a aAt 
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( 25 ) 


where < > means the average of the quantity within the brackets in 
the direction of the subscript on the last bracket. The stretched 
coordinate and the salt continuity equation could have been shown 
also but these add nothing to the discussion at this point. If 
h = zmax, a » Az, these equations are identical to the Interior 
equations. If the surface does not coincide with a z grid point 
(as it never would), the uJh" 0* is the normal behavior of 

these equations. For convenience, a certain number of z-grid points, 
zmax, are always used, in which case a might be negative. This 
feature imposed the condition that "a" never excedes ] 

in magnitude. This creates no problem; it is merely described so 
that the logic of the calculation may be understood. 


Pressure, the last dependent variable to be calculated in the 
sequence of each time step, can be computed by a spatial marching 
procedure similar to that used to calculate U* . At the free surface 
the pressure will be known, and pressure at loser grid rows can be 
computed from an integrated form of the z-momentum equation. An 
initial adjustment must be made to the pressure at the top grid 
plane to account for the free surface height. 


Boundary Conditions 


Boundary conditions for each dependent variable must be enforced 
after the sequence in which the variable is computed during the time 
step. Boundary conditions were applied as shown in Figure 6. 

A no- slip condition was enforced along the bottom by forcing 
the velocity at grid points adjacent to the bottom to fit a boundary- 
layer-type profile. Grid points lying under the bottom were by- 
passed. The vertical shear component was set equal to zero as a 
boundary condition. Wind shear at the surface was empirically 
related to wind velocity. Heat flux to and from the atmosphere, 
although neglected in this analysis, could have been Included as 
surface boundary conditions. 


Reflection principles were used along the bay banks, where it 
was felt that the coarse horizontal grid spacing was insufficient 
to enforce a no-slip condition. 

Velocities through the tidal inlets were obtained by extrapo- 
lation, for a given tide level and cut-point on the salinity dis- 
tribution. The cut-point is defined as a level below which there 
is pure salt water. These conditions could be held constant with 
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time to obtain a steady-state solutiont or varied as a known func- 
tion of time to obtain a time-dependent solution. The tidal effects 
will be modeled by Imposing time varying boundary condition’* at the 
seaward inlets to the harbor. Over part of the tide cycle tne river 
flows will also have to be time dependent. Flows to and from marshes 
should also be Included as part of the specified boundary conditions. 
These latter two effocts have not yet been included in the model. 


RESULTS 


The equations presented in the previous chapter cons>"l tute the 
mathematical model of Mobile Bay, The solution technique , •'•rented 
has been implemented with a computer program; this progra- la listed 
in the Appendix. Mathematically, the model is a set of partial dif- 
ferential equations which are solved for a particular set of auxiliary 
conditions. Since the desired solution is a function of time, ini- 
tial conditions throughout the bay are required; these are not 
available. In fact, one of the major benefits of using a model is 
to first predict a quasi-steady or a periodic motion such that an 
accurate statement of Initial conditions is not required. A pro- 
posed scheme of generating the necessary initial conditions is 
presented. The speed and economy of obtaining a solution with this 
method can now be determined. Initial calculations indicate that 
more Improvement in efficiency is warranted before extensive cal- 
culations are made. Several methods of accomplishing this improve- 
ment are outlined. 

The concept for modeling Mobile Bay is that the three-dimensional, 
unsteady conservation laws are solved for a set of boundary condi- 
tions and an arbitrary set of assumed initial conditions. The 
boundary conditions are held fixed such that the solution to the 
equation is driven to an asymptotic limit. This limit then becomes 
an Initial value for a truly unsteady solution. This idea was pre- 
sented in (8) and applied to modeling a bayou flow problem in (9). 

The river flows are held fixed. The tide level and a cut-point 
below which there is only sea water is fixed at the Mobile Point 
entrance to the bay. Flow through the Cedar Point - Dauphin Island 
entrance is obtained by extrapolation. These boundary conditions 
and an initial estimate of all dependent variables are sufficient 
to begin the sequence of calculations. 

The program to perform these calculations is listed as CHANNEL. 
CHANNEL describes the bay for the stated type of boundary conditions 
and, in addition, places a specified wind shear on the surface and 
a no slip condition on the bottom. CHANNEL is a combination of MOBILE 
2 and the bayou program from ref. (9). The bayou calculation is 
Included as a description of the ship channel. MOBILE 2 * is also 

*The 2 in its name means it is the second grid configuration which 
has been used. Generally, there is no difficulty jn changing the 
number of grid points in a calculation, but for thi.s bay calcula- 
tion some care is necessary, because internal program logic requires 
that adjacent boundary grids differ by no more than one. This means 
that the ratios of grid points used should be kept roughly the same. 
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a three-dimensional, unsteady bay model, but the somewhat more sim- 
ple bottom location specification allows an overall more simple 
solution. This development of programs by addition of pre-tested 
subroutines to a basic main program Is by design; this allows one 
to test and Incorporate '* features Into the calculation without 
undue program rewriting. 

It should be noted that CHANNEL is the culmination of several 
previous studies which has resulted in significant Improvements in 
this basic analytical approach. Specifically, these improvements 
are: 


1) An equation set which has been optimized to run faster by 
eliminating terms which an order of magnitude analysis showed to 
be of negllble magnitude. 

2) A windward differencing scheme was used since it appears 
to be the fastest explicit scheme yet devised. 

3) Wind-shear is included as a surface boundary condition. 

4) The method of applying tidal conditions has been proved 
to be adequate on the more simple bayou analysis. 

Starting with these program features a grid system was devel- 
oped for Mobile Bay as shown in Fig. 6. Calculations were then 
initiated to describe the Bay. 

A zero mean- sea- level tide condition and a set-point corres- 
ponding to field observations at Mobile Point were used as bound- 
ary conditions to seek an asymptotic solution. The constancy and 
smallness of the program generated error functions indicate that 
suitable progress has been made to begin seeking an unsteady solu- 
tion. The error functions (SEAR, UBAR, VBAR) are gross measures 
of the changes in S, U,. Uy between computation steps. 

Study of the calculations and review of som^ of the development 
steps necessary to obtain these calculations yield the following 
results. 

Stable calculation time is about 0.5 seconds. Modest changes 
in the grid system point density and using MOBILE 2 rather than 
CHANNEL does not appreciably change this step size; however, both 
of these changes do affect the total amount of time necessary to 
perform a given number of calculation steps. The significant factor 
which must be considered is the amount of real time which is modeled 
in a unit of computation time. For CHANNEL and MOBILE 2, these 
numbers are, about 6~^ , for grid system number 2 (12 x 20 x 9). 
(Actually, the nine is an eight in MOBILE 2.) Grid system 1 (17 x 
30 X 15) gives about 25"^ as the real time to computing time ratio. 

Data at intermediate numbers of calculation steps reveal that 
velocity and surface elevation fields establish before concentration 
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fields. Afi expected for the case run the velocities are small, 
probably because the major driving force in creating flow is the 
tidal oscillation which have not yet been calculated. 

Initially guessed fields take a long time to change - of the 
order of hours of real time. 

The most significant result is that the program has been 
developed, is running without apparent errors, and is responding 
as expected to applied boundary conditions. 

An unsteady solution has not yet been sought, mainly because 
the program has only now been developed and run to obtain these 
initial conditions. Secondly, several observations can be made on 
the calculations and available experimental data which have a bear- 
ing on the Immediate steps which need to be made in researching this 
problem. 

River flows as boundary conditions need to be specified so that 
they can reverse. This means specifying them as a function of time 
during a tide cycle. Experimental field data at one river stage 
condition are available and can be used for such a specification. 

The ratio of real time to computer time needs to be improved 
by some or all of the following means. 

1) Determining if using a larger number of time steps in the 
species equation than in the momentum equations Improves computational 
efficiency. 

2) Running on a faster computer. Some runs on machines avail- 
able to NASA would expedite this study. 

3) Transform the time variable to reduce the real time length 
of a tide cycle. Reductions in cycle from 24 hours to about 4 (the 
lag time for a tidal flucturation to travel from the south to the 
north end of the Bay) should prove feasible. 

4) Devise an improved computation procedure. Possibilities 
from the simple expedient of devising a new guessing procedure for 
initialization to using a new differencing formulation will be 
considered. 

In summary, a mathematical model of Mobile Bay has been devel- 
oped. This model is now ready to use for parametric studies to 
predict the bay's dynamic behavior. Economy dictates that the exact 
procedure for accomplishing the required computation be carefully 
planned and executed. 
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APPENDIX 


plBcusalon of the Program ; This appendix deacribea a computer pro- 
gram which can be used Sot modeling a bay or eatuaiy with a shipping 
channel. The mathematical model used In this program Is a three- 
dimensional, time-dependent flow field model which can be used to 
describe the type flow found In any bay or estuary by changing 
appropriate boundary conditions and geometry. 

Mobile Bay was selected as the bay to be modeled. The compu- 
tational techniques and equations solved are presented In the 
chapter on model techniques of the main paper. Comment cards are 
used throughout In the computer program to Indicate which operations 
are being performed, A simplified flow sheet Is presented In 
Figure 1. 

Input Guide ; All data cards to be used in the program are read by 
SUBROUTINE PRELIM, Proceeding eaciA READ statement are comment cards 
which define each term to be read and specify the appropriate units. 
These data are nondlmenslonallzed before using them for computations. 
Two basic formats are used for Input, IlO and FlO.5. The IlO Is an 
Integer format consisting of 10 columns right justified; the FlO.5 
Is a floating point format occupying 10 columns with a decimal point 
punched on the card. 

At the completion of the run, data are written on a tape in 
SUBROUTINE ECRIRE. This provides data for a restart capability or 
a plot program. 

When the restart capability Is used, data are read from a tape 
by SUBROUTINE LIRE. The option of calling this subroutine Is con- 
trolled by a data card as described In SUBROUTINE PRELIM. 

Table 1 provides the details of the card Input. 

Description of Printout ; This section presents a description of 
the CHANNEL program output. 

Initially, the Input parameters are printed in order to identify 
the problem and to verify that correct values were used. Then a 
complete description of the grid systems In the three coordinate 
directions are printed. This description Includes index numbers and 
dimensionless and dimensional values for the locations of the grid 
planes as well as stretching function values where applicable. 

The physical location of the boundaries of the bay are printed 
as It was read Into the program and as the program describes it In 
terms of the above grid systems. If desired, a description of the 
floor of the bay can be printed. This portion of the printout is 
generally suppressed by branching around it with a G0 T0 statement 
because of Its excessive length and low usage. This same comment 
also applies to a portion of the printout that describee the heights 
of the banks along the boundaries of the bay. Then the maximum 
velocities, east and west bank locations and depths of the rivers 
entering the north end of the bay are printed. 
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Next, terms monitoring the overall changes In the salinities 
and horizontal components o£ velocity £or a single time step are 
printed £or selected time steps. If the number o£ time steps sped” 
fled In the Input data Is achieved, a message stating that a stable 
run was made Is given and then values o£ the flow field variables 
are printed at selected Y-Z grid planes^ with velocities in M/S, 
salinities as fractions of Gulf water salinity and pressures In non" 
dimensional form. If the surface becomes too rough or a negative 
pressure below the surface Is encountered, a message Identifying 
the cause for premature termination along with Its location In the 
grid system is given. Values of the flow field variables are then 
printed as above. Finally, the free surface heights relative to 
the top horizontal grid plane are printed In CM in the form of a 
map for ease of Interpretation. 

Program Listing ; The following pages contain a listing of CHANNEL 
and the required input data. 


1 Which Y-Z grid planes are selected is controlled by LOGICAL IF 
statements that can be modified as needed. 
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ncuiE 1 


SIMPLinED FLOW CHART FOR CHANNEL 


C start) 

\READ INPUT PARAMETERS/ 

\ f 

PERFORM PRELIMINARY CALCUUTIONst 

1 

\PRINT INPUT PARAMETERS/ 


YES 



NO 


START. 


INITIALIZE FLOW 
FIELD VARIABLES 


READ FLOW FIELD 
^VARIABLES FROM TAPE 







7 


SOLVL X MON. EQN. FOR U AND 
SALT CONS. EQN. FOR S IN CHANNEL 
AT NEW TIME 

* • ~ ~ 
SOLVE VOLUME DILATATION 
FOR W IN CHANNEL AT NEW TIME 
* 

I SOLVE X & Y MOM. E^IS. FOR U & V 
AMD SALT CGNS. EQN. FOR S 
THROUGHOUT BAY AT NEW TIME 

« 

SOLVE VOLUME DIIATATION E(|l. 
FOR W THROUGHOUT BAY AND 
FOR SURFACE HEIGHT AT FREE 
INTERFACE AT NEW TIME 


\ 

NO 

1 SOLVE Z MOM, 

, EQN. FOR P 1 

1 THROUGHOUT BAY AT NEW TIME | 



WRITE FLOW FIELD 

Wariables on tape 


SOLVE Z MOM. EQM. FOR Pl 
IN CHANNEL AT NEW TIME 



I 


PRINT FLOW 
lELD VARIABLES 

I 


7 


(stop) 
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TA.BLE 1 

INPUT PARAMETERS FOR CHANNEL 


card no. 

VARIABLES 

FORMAT 

1 

LTAPE 

110 

2 

IMAX, JMAX, KMAX, JWAG, KMXC 

110 

3 

NMAX, NPRINT, NEAR 

110 

4 

FETCH, YMAX 

F10.5 

5 

XBANKN, XBANKS 

F10.5 

6 

DTREAL, YK2, ACCEL 

F10.5 

7 

HREF, SIGMAT, PERI0D, TIDE 

Flo. 5 

8 

ROUGH, V0NKAR 

FlO.5 

9 

VWIND, DIRECT, PHI 

Flo. 5 

10 

CDEPTH, CUTPT 

F10.5 

11 

UMAX 

110 

12* 

S0UTH(II), EAST(II), WEST(II) 

FlO.5 

13 

C0NVRT 

F10.5 

14** 

YW(I), YE(1), DEPTH(I), URIV(I), RIV 

FlO.5 


* UMAX of these cards are read. 

** A card is read for each of the 4 rivers entering the north end 
of the bay. 
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5JMl=S(MO« I « JMl ,KBOT) 
SJPl=S(MO, 1 ,JP1 ♦KBOT> 
VJMl = VCHO, I , JMl .K0CT) 
VJPl = V(MO. I * JP 1 tKBOT> 
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S(MN« I, J«K) = SNEW 

SBAR = SBAR+ABS( SNEK»-Sai.O> 
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JE=JEAST( IMAX) 

DO 34 J=JW.JE 
KBOT=KFLOOR( IMAX« J) 
DO 34 K=KBOT«KMAX 


S<MN* I*J.K)=S(MNtIMAXM1 . J.K) 
34 CONTINUE 

DO 36 I=2»IMAX 
J«=JHESTf 1 » 

JE=JEAST( I } 
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SUBROUTINE UVSNEV* 

THIS SUBROUTINE CALC. NE»* VALUES OF U.V.C S 
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30 CONTINUE 
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40 CONTINUE 
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COMMON/FLOOR/KFLOORC 12.20) *ZB(12.20) .KFLORCC 12) .2HC( 12) 
COHMUN/PULL./X( 12)«Y(20)«Z(ll) . SY (20 ) . SYY ( 20 ) .YK2.2C( 1 1 ) 
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1 853 


!•! 1*35 e«53 0.213 13.0 

4.4 4.63 6.10 0.166 17.0 

5.5 5.6 7.92 0.229 16.0 
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